#include <string>
#include <algorithm>
#include <iostream>
#include <vector>
#include <sstream>
#include <boost/lexical_cast.hpp>
#include "components/portable_timer.h"
#include "components/progress_indicator.h"
#include "components/pg_db_cursor.h"
#include "get_arg.h"
#include "components/bl2seq_mine.h"
#include "components/perlutil.h"

using std::cerr;
using std::cout;
using std::string;
using std::vector;
using boost::lexical_cast;

struct t_gi_seq
{
	string gene_id;
	string acc;
	string pi;
	string seq;
	t_gi_seq(const string& g, const string& a, const string& p, const string& s):
				gene_id(g), acc(a), pi(p), seq(s){}
};

typedef vector<t_gi_seq> t_vec_gi_seq;

//________________________________________________________________________________________

//		get_species_list

//________________________________________________________________________________________
bool get_species_list(t_pg_cursor& pg_db, vector<string>& vec_species)
{
	string SQL_command = "SELECT taxon FROM genomes.ensembl_parameters";
	if (!pg_db.ExecTuplesOk(SQL_command))
	{
		pg_print_sql_err("retrieve the list of species from postgreSQL",
							pg_db, SQL_command, cerr);
		return false;
	}
	unsigned cnt_species = pg_db.Tuples();
	while (cnt_species)
		vec_species.push_back(pg_db.getvalue(--cnt_species, 0));

	return true;
}



//________________________________________________________________________________________

//		get_gene_clusters

//________________________________________________________________________________________
bool get_gene_clusters(t_pg_cursor& pg_db, const vector<string>& species,
						const unsigned cluster_id,
						t_vec_gi_seq& vec_gi_seq,
						std::ostream& VERBOSE)
{

	string cluster_id_str = lexical_cast<string>(cluster_id);
	string sql_cmd;
	for (unsigned i = 0; i < species.size(); ++i)
	{
		if (i)
			sql_cmd += " UNION ";
		sql_cmd += "SELECT gene_id, acc_code, pi, sequence "
					"FROM orthologs.tribe_clusters "
					"NATURAL JOIN " + species[i] + ".ens_id "
					"NATURAL JOIN sequences "
					"WHERE cluster_id = " + cluster_id_str;

	}
	sql_cmd += " ORDER BY gene_id ";
	{
		// create binary cursor
		if (!pg_db.Declare(sql_cmd, false))
			return pg_print_sql_err("create a cursor to "
									"retrieve sequences for cluster " + cluster_id_str,
									pg_db,
									sql_cmd, cerr);

		unsigned batch_size = 50;
		// get data in batches
		while (pg_db.Fetch(batch_size))
		{
			int curr_batch_size = pg_db.Tuples();
			if (!curr_batch_size)
				break;
			for (int i = 0; i < curr_batch_size; ++i)
				vec_gi_seq.push_back(t_gi_seq(	pg_db.getvalue(i, 0),
												pg_db.getvalue(i, 1),
												pg_db.getvalue(i, 2),
												pg_db.getvalue(i, 3)));

		}

		// close cursor
		if (!pg_db.Close())
			return pg_print_sql_err("close the  cursor for "
									"retrieving sequences for cluster " + cluster_id_str,
									pg_db, "CLOSE CURSOR",
									cerr);
	}
	return true;
}


/*
bool expand_gaps(const string& gaps_str, string& prot)
{
	// empty string encodes no gaps!
	if (! gaps_str.length())
		return true;


	// split string of gaps
	vector<string> gaps;
	perl_split(gaps_str,  "|",  gaps);
	for (int i = gaps.size()-1; i >= 0; --i)
	{
		string::size_type pos = gaps[i].find(':');
		if (pos == string::npos)
			return false;
		// gaps should not come after a protein or be humungously long
		unsigned beg = atoi(gaps[i].c_str());
		unsigned len = atoi(gaps[i].c_str() + pos + 1);
		if (beg > prot.length() || len > 80000)
			return false;
		prot.insert(beg, string(len,  '-'));
	}
	return true;
}
*/


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_portable_timer timer;

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// saving parameters
	string exe_parameters;
	for (int i = 0; i < argc; ++i)
		exe_parameters += string(argv[i]) + ' ';

	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	// get command line arguments
	t_program_args args;
	if (!GetArguments(argc, argv, args))
		return 1;
	assert(args.db_origin.length() > 0 && args.db_origin.length() <= 5);


	//____________________________________________________________________________________
	// 		Connect to PostgreSQL
	//____________________________________________________________________________________
	t_pg_cursor pg_db(args.db_login.c_str(), "panda_cursor");	//"dbname=test user=postgres"
	if ( pg_db.ConnectionBad() )
	{
		cerr	<< "\nERROR!! Connection to database ["
							<< args.db_login << "] failed.\n"
							<< pg_db.ErrorMessage() << "\n\n";
		return 1;
	}

	//____________________________________________________________________________________
	// 		get largest cluster id
	//____________________________________________________________________________________
	string SQL_command = "SELECT MAX(cluster_id) FROM orthologs.tribe_clusters";
	if (!pg_db.ExecTuplesOk(SQL_command))
	{
		pg_print_sql_err("retrieve the maximum tribe cluster id from postgreSQL",
							pg_db, SQL_command, cerr);
		return 1;
	}
	unsigned max_cluster_id = lexical_cast<unsigned>(pg_db.GetValue(0, 0));
	args.last_cluster = std::min(args.last_cluster, max_cluster_id+1);
	args.first_cluster = std::min(args.first_cluster, max_cluster_id);



	//____________________________________________________________________________________
	// 		Get ENSG to species mapping
	//____________________________________________________________________________________
	args.VERBOSE() << "Getting species list\n";
	std::vector<string> vec_species;
	if (!get_species_list(pg_db, vec_species))
		return 1;

	// prepare blast
	//BLAST_OptionsBlkPtr options = create_bl2seq_options();
	t_blast_options options(false, 1e-5, "BLOSUM90");

	//____________________________________________________________________________________
	// 		Blast all against all for each cluster
	//____________________________________________________________________________________
	args.VERBOSE()	<< "Do all against all blast for clusters "
					<< args.first_cluster	<< " to "
					<< args.last_cluster	<< "...\n";
	unsigned cnt_bl2_seq_ops = 0;
	{
		t_progress_indicator dots(args.VERBOSE(), 50);
		for (unsigned clust_id = args.first_cluster; clust_id < args.last_cluster; ++clust_id)
		{
			t_vec_gi_seq vec_gi_seq;
			if (!get_gene_clusters(pg_db, vec_species, clust_id, vec_gi_seq, args.VERBOSE()))
				return 1;
			for (unsigned i = 0; i + 1 < vec_gi_seq.size(); ++i)
				for (unsigned j = i + 1; j < vec_gi_seq.size(); ++j)
				{
					if (vec_gi_seq[i].gene_id == vec_gi_seq[j].gene_id)
						continue;
					++dots;
					double			e__value,	bitscore;
					vector<t_gap>	gaps1, 		gaps2;
					unsigned 		align_beg,	align_len, total_len;
					if (!do_bl2seq(	options,
									vec_gi_seq[i].seq,
									vec_gi_seq[j].seq,
									e__value,	bitscore,
									gaps1,	gaps2,
									align_beg,	align_len,
									total_len))
					{
						args.CERR() << clust_id << ":"
									<< vec_gi_seq[i].gene_id	<< "\t"
									<< vec_gi_seq[j].gene_id    << "\t"
									<< vec_gi_seq[i].acc        << "\t"
									<< vec_gi_seq[j].acc        << "\t"
									<< vec_gi_seq[i].pi         << "\t"
									<< vec_gi_seq[j].pi         << "\n"
									<< vec_gi_seq[i].seq        << "\n"
									<< vec_gi_seq[j].seq        << "\n";
						continue;
					}
					string		gap_str1 = gap_array_to_str(gaps1);
					string		gap_str2 = gap_array_to_str(gaps2);
					args.f_output_file
										<< clust_id					<< "\t"
										<< vec_gi_seq[i].gene_id	<< "\t"
										<< vec_gi_seq[j].gene_id    << "\t"
										<< vec_gi_seq[i].acc        << "\t"
										<< vec_gi_seq[j].acc        << "\t"
										<< vec_gi_seq[i].pi         << "\t"
										<< vec_gi_seq[j].pi         << "\t"
										<< e__value				    << "\t"
										<< bitscore         	    << "\t"
										<< align_beg       	 	    << "\t"
										<< align_len       	 	    << "\t"
										<< total_len       	 	    << "\t"
										<< (total_len ? (align_len * 1.0 / total_len) : 0.0)
																	<< "\t"
										<< gap_str1         	    << "\t"
										<< gap_str2         	    << "\n";
					string seq1 = vec_gi_seq[i].seq;
					string seq2 = vec_gi_seq[j].seq;
					expand_gaps_prot(gaps1, seq1);
					expand_gaps_prot(gaps2, seq2);
					args.f_output_file << seq1 << "\n" << seq2 << "\n";
					args.f_output_file
									<< string(seq1, align_beg, align_len) << "\n"
									<< string(seq2, align_beg, align_len) << std::endl;
					cnt_bl2_seq_ops++;
				}


		}
	}

	//free_bl2seq_options(options);
	cerr << timer.reset() << " Seconds for " << cnt_bl2_seq_ops << " iterations\n";

	return 0;
}

